import json
import math
import os
import requests
import cv2
import numpy as np
from osgeo import gdal, osr, ogr
import shutil
import time
import sys


def initSetting():
    global setting
    setting = getConfig('setting.json')
    mkdir('.\\'+setting['tilePath'])
    mkdir('.\\'+setting['outputPath'])


def swap(a, b):
    a, b = b, a
    return a, b

class Point:
    def __init__(self, lon, lat):
        self.lon = lon
        self.lat = lat

class Box:
    def __init__(self, point_lt, point_rb):
        self.point_lt = point_lt
        self.point_rb = point_rb


def build_url(x, y, z):
    return setting['url'].format(x=x, y=y, z=z)

def download(x, y, z, path):
    proxies = setting["proxies"]
    url = build_url(x, y, z)
    path = path + "\\{z}\\{x}\\".format(z=z, x=x)
    if not os.path.exists(path):
        os.makedirs(path)
    filepath = path + "{y}.tif".format(y=y)
    if os.path.exists(filepath) and os.path.getsize(filepath) > 400:
        pass
    else:
        for x in range(0,3):
            # response = requests.get(url, proxies=proxies)
            response = requests.get(url)
            if response.status_code == 200:
                with open(filepath, "wb") as f:
                    f.write(response.content)
                break;
            else:
                print("network error!")

def xyz2lonlat(x, y, z):
    n = math.pow(2, z)
    lon = x / n * 360.0 - 180.0
    lat = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    lat = lat * 180.0 / math.pi
    return lon, lat

def lonlat2xyz(lon, lat, zoom):
    n = math.pow(2, zoom)
    x = ((lon + 180) / 360) * n
    y = (1 - (math.log(math.tan(math.radians(lat)) + (1 / math.cos(math.radians(lat)))) / math.pi)) / 2 * n
    return int(x), int(y)

def cal_tiff_box(x1, y1, x2, y2, z):
    LT = xyz2lonlat(x1, y1, z)
    RB = xyz2lonlat(x2 + 1, y2 + 1, z)
    return Point(LT[0], LT[1]), Point(RB[0], RB[1])

def core(path,lt,rt,z):
    point_lt = Point(lt[0], lt[1])
    point_rb = Point(rt[0], rt[1])
    x1, y1 = lonlat2xyz(point_lt.lon, point_lt.lat, z)
    x2, y2 = lonlat2xyz(point_rb.lon, point_rb.lat, z)
    count = 0
    all = (x2-x1+1) * (y2-y1+1)
    for i in range(x1, x2+1):
        for j in range(y1, y2+1):
            if(count%100==0):
                print("{m}/{n}".format(m=count, n=all))
                print("\r", end="")
                print("进度: {}%: ".format(int(count/all *100 )), end="")
                os.system('cls')
            download(i, j, z, path)
            count += 1
    print('该幅图下载完毕进入合并')
    merge(x1, y1, x2, y2, z, path)

def merge(x1, y1, x2, y2, z, path):
    row_list = list()
    for i in range(x1, x2+1):
        col_list = list()
        for j in range(y1, y2+1):
            col_list.append(cv2.imread(path + "\\{z}\\{i}\\{j}.tif".format(i=i, j=j,z=z)))
        k = np.vstack(col_list)
        row_list.append(k)
    result = np.hstack(row_list)
    cv2.imwrite(path + "\\merge.tif", result)
    print('合并完成')


def registration(input_path, out_path, top_left, bottom_right, ik, jk, name,srs):
    """
    基于python GDAL配准
    :param input_path: 需要配准的栅格文件
    :param out_path: 输出配准后的栅格文件位置
    :param top_left: 左上角坐标
    :param bottom_right: 右下角坐标
    :param ik: 行空白分辨率
    :param jk: 列空白分辨率
    :return:
    """
    # 读文件
    YG2_dataset = gdal.Open(input_path)  # 打开文件，用这个tif图的数值信息
    x, y = YG2_dataset.RasterXSize,YG2_dataset.RasterYSize  # 栅格矩阵的列数
    del YG2_dataset
    # 打开栅格文件
    dataset = gdal.Open(input_path, gdal.GA_Update)
    # 构造控制点列表 gcps_list
    gcps_list = [gdal.GCP(top_left[0], top_left[1], 0, 0, 0),
                 gdal.GCP(bottom_right[0], top_left[1], 0, x - jk, 0),
                 gdal.GCP(top_left[0], bottom_right[1], 0, 0, y - ik),
                 gdal.GCP(bottom_right[0], bottom_right[1], 0, x - jk, y - ik)]
    # 设置空间参考
    spatial_reference = osr.SpatialReference()
    if srs == 4528:
        spatial_reference.SetWellKnownGeogCS('CGCS2000')
    else:
        spatial_reference.ImportFromEPSG(srs)
    # 添加控制点
    dataset.SetGCPs(gcps_list, spatial_reference.ExportToWkt())
    # tps校正 重采样:最邻近法
    layer = 'memory\\'+name+".shp"
    dst_ds = gdal.Warp(out_path, dataset, format='GTiff', tps=True, width=x, height=y,
                       resampleAlg=gdal.GRIORA_NearestNeighbour,cutlineDSName=layer)
    del dataset

def mkdir(path):
    folder = os.path.exists(path)
    if not folder:  # 判断是否存在文件夹如果不存在则创建为文件夹
        os.makedirs(path)  # makedirs 创建文件时如果路径不存在会创建这个路径
def getConfig(filePath):
    with open(filePath,"r",encoding="UTF-8") as jsonFile:
        config = json.loads(jsonFile.read())
    return config


def ReadShp(_filename):
    # 支持中文路径
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
    # 支持中文编码
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    # 注册所有的驱动
    ogr.RegisterAll()
    # 打开数据
    ds = ogr.Open(_filename, 0)
    if ds == None:
        return ("打开文件失败！")
    # 获取数据源中的图层个数，shp数据图层只有一个，gdb、dxf会有多个
    # iLayerCount = ds.GetLayerCount()
    # 获取第一个图层
    oLayer = ds.GetLayerByIndex(0)
    if oLayer == None:
        return ("获取图层失败！")
    # 对图层进行初始化
    oLayer.ResetReading()
    # 输出图层中的要素个数
    num = oLayer.GetFeatureCount(0)
    layer_list = []
    # 获取要素
    sr = osr.SpatialReference()
    sr.ImportFromEPSG(setting['src'])
    for i in range(0, num):
        ofeature = oLayer.GetFeature(i)
        fieldName = ofeature.GetFieldAsString('NAME')

        driver = ogr.GetDriverByName('ESRI Shapefile')
        data_source = driver.CreateDataSource('memory')
        layer = data_source.CreateLayer(fieldName,sr)
        layer.CreateFeature(ofeature)

        bbox = layer.GetExtent()
        layer_list.append([fieldName,[bbox[0],bbox[3],bbox[1],bbox[2]]])
    ds.Destroy()
    del ds
    return layer_list

if __name__ == '__main__':
    if os.path.isdir('memory'):shutil.rmtree('memory')
    initSetting()
    # 存在多个版本投影相关proj文件，需要进行指定
    if(setting['proj']):os.environ['PROJ_LIB'] = setting['proj']
    # # 设置名称和范围
    info = ReadShp(setting['shapefile'])
    n = 1
    time_start = time.time()  # 记录开始时间
    for z in range(setting['startlevel'],setting['endlevel']+1):
        mkdir('.\\'+setting['outputPath']+'\\' + str(z))
        for i in info:
            # 影像下载
            print('当前层级：{z}当前进度:{n}，总影像数目：{lens}'.format(z=z,n=n, lens=len(info)))
            path = '.\\'+setting['tilePath']+'\\' + i[0]
            mkdir(path)
            core(path, (i[1][0], i[1][1]), (i[1][2], i[1][3]), z=z)  # 调整下载级别
            registration(path + '\\merge.tif', '.\\'+setting['outputPath']+'\\'+str(z)+'\\' + i[0] + '.tif', (i[1][0], i[1][1]),
                         (i[1][2], i[1][3]), setting['revolution'], setting['revolution'], i[0], setting['src'])
            time_end = time.time()
            n+=1
        n = 1
    print("共花费{s}秒".format(s=str(time_end - time_start)))
    if os.path.isdir('memory'):shutil.rmtree('memory')




    # http://khms0.google.com/kh/v=893?&x={3399}&y={1682}&z={12}